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Abstract 

We study compressible fluid flow in narrow two-dimensional channels using a novel 
molecular dynamics simulation method. In the simulation area, an upstream source 
is maintained at constant density and temperature while a downstream reservoir is 
kept at vacuum. The channel is sufficiently long in the direction of the flow that the 
finite length has little effect on the properties of the fluid in the central region. The 
simulated system is represented by an efficient data structure, whose internal elements 
are created and manipulated dynamically in a layered fashion. Consequently the code 
is highly efficient and manifests completely linear performance in simulations of large 
systems. We obtain the steady-state velocity, temperature, and density distributions 
in the system. The velocity distribution across the channel is very nearly a quadratic 
function of the distance from the center of the channel and reveals velocity slip at the 



boundaries; the temperature distribution is only approximately a quartic function 
of this distance from the center to the channel. The density distribution across the 
channel is non-uniform. We attribute this non-uniformity to the relatively high Mach 
number, approximately 0.5, in the fluid flow. An equation for the density distribution 
based on simple compressibility arguments is proposed; its predictions agree well with 
the simulation results. Validity of the concept of local dynamic temperature and the 
variation of the temperature along the channel are discussed. 

PACS numbers: 47.40.Dc and 47.60.+i 
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I. Introduction 



The technique of molecular dynamics (MD) simulation has been widely used to 
study non-equilibrium fluids. Because of limitations imposed by finite computational 
capacity with regard to both memory and speed, this method has been used princi- 
pally to determine the behaviors of fluid systems on time and distance scales within 
a few orders of magnitude of r c and r c , respectively. Here, r c ~ 1CP 13 sec is the 
collision duration and r c ~ 1CP 8 cm is a molecular size.0 Only molecular properties 
of the fluid can be obtained in this range of time and distance, and that is also the 
domain of experimental neutron scattering measurements. Realistic examination of 
the hydrodynamic properties of flow is still beyond the reach of most molecular dy- 
namics simulations, although there are many attempts to simulate larger systems on 
longer time scales. Among these are the work of Koplik, et al., 2,3 Hannon, et al.$ 
and Bhattacharya and Lie.!'! In particular, Hannon, et al., have obtained velocity 
and temperature distributions across channels in which flow is occurring. Their re- 
sults agree well with simple hydrodynamic predictions for incompressible fluid flow. 
Also, Bhattacharya and Lie have obtained similar velocity profiles and have computed 
boundary slip coefficients. 

These simulations of fluid flow have a common property. Periodic boundary con- 
ditions are introduced along the flow direction in the interest of reducing the amount 
of computation needed to obtain useful results. The imposition of translational in- 
variance along the flow direction makes it necessary to introduce a "gravitational" 
field to induce flow. In order to induce appreciable flow, this field must be given a 
strength much larger than the earth's field g, for example, as large asi 10 12 g. As a 
consequence of the gravitational field, regular rescaling of the particles' kinetic en- 
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ergies is required; hence one cannot reliably study properties of the fluid having to 
do with energy or heat flow. Also, it is not possible to study variations of the flow 
properties along the flow direction. 

In this paper, we present a new method of simulation in which the channel is 
of finite length in the flow direction without periodic boundary conditions. At one 
end of the channel is a source region which is maintained at constant density and 
temperature by introducing particles as needed. At the other end is a sink region 
which is maintained at vacuum; that is, any particle which moves into this region is 
removed from the system. Hence the pressure or density gradient along the channel 
is primarily responsible for instigation of the flow. We also use a novel layered data 
structure to improve the efficiency of the code and the utilization of storage. Our 
methods make it possible to simulate systems containing 20,000 or more particles for 
more than 10 6 time steps on a DECstation 3100. 

Section II contains a description of the system simulated and of our numerical 
techniques. The results are in Sec. Ill, and Sec. IV contains a discussion and conclu- 
sions. 

II. Model and numerical methods 

The geometry of the system is indicated in Fig. 1. The channel has a length L in 
the x direction and width w in the y direction. Flow is along the x direction with the 
source region at < x < L\ and the vacuum or sink located at x > L. The channel is 
closed at x = so that particles cannot escape into the region x < 0. Typical channel 
sizes that we have investigated are L = 400a, w = 100cr, and L\ = 100cr, where a is 
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the particle size parameter (diameter) in the Lennard-Jones potential 



V(r) = 4e 



a\ 12 {cj^ 



(1) 



which we have employed for the interparticle interactions. The interaction is trun- 
cated at r > 2<j .i In addition to e and a, the only other parameter describing the 
properties of the molecules is the mass m. We have used values appropriate for ar- 
gon, i.e., a = 3.4rA, e/k = 119.76 K, and m = 6.67 x 10~ 23 g; k is the Boltzmann 
constant. We also introduce a basic time constant r = ^/mcr 2 /48e ~3x 10 -13 sec. 

Our procedure for handling collisions of the particles with the walls, which are 
the surfaces x = 0, y = 0, and y = w, is to give the recoiling particles a Maxwell- 
Boltzmann velocity distribution on a half-space. For example, after colliding with 
the wall at x = a particle is given a velocity v = v x x + v y y (where v x > 0) 
with a probability proportional to exp[— m(v x + Vy)/2kT] where T is the specified 
wall temperature. The molecular dynamics simulation proceeds in the conventional 
manner except that when the particle density in the source region drops by a small 
amount, typically after several hundred time steps, we inject particles into this region 
in order to bring the density here back up to some prescribed value. That is done 
by, first, letting the end wall (x = 0) act as a piston and uniformly compress the 
particles in the source region (but not the remainder of the channel) by a sufficient 
amount to bring the density in this compressed volume up to the desired level, and, 
second, restoring the end wall to its original position and injecting an appropriate 
number of particles into the empty space next to the wall. These added particles 
are given a Maxwell velocity distribution at the same temperature T as that of the 
walls. By making this adjustment to the number of particles sufficiently frequently, 
we only have to inject a few particles each time and the relative amount by which the 



source region is compressed is very small. By comparison, the entire system typically 
contains ten to twenty thousand particles, and the source region is one- fourth of the 
channel. 

The initial state of the system is always taken to have some preset density of 
particles with a Maxwell-Boltzmann velocity distribution and no additional particles 
elsewhere. Thus, as time proceeds, particles make their way down the channel, driven 
by the pressure or density gradient, and are removed at the far end. After a sufficiently 
long time a steady state is achieved and various quantities such as the distributions 
of flow velocity, particle density, particle current density, and temperature, may be 
computed. We typically run the simulation in the steady state for some 10 6 time 
steps with one time step being St = 3 x 10~ 2 t ~ 9 x 1CP 15 sec. In addition, some 10 6 
time steps are needed for the system to reach the steady state. 

The most critical part of a molecular dynamics simulation is the force calculation! 
which would require 0(N 2 ) steps in a brute force calculation if N is the number of 
particles. Given short-ranged interactions, such as the Lennard-Jones potential, one 
typically uses Ver let's neighbor listJH to lower the number of steps to 0(NN c ), where 
N c is the maximum number of particles that can fit within the range of a particle's 
potential. This method has high memory requirements and also requires a number 
of steps 0(N 2 ) to update the neighbor list. A superior method! in large systems is 
to divide the simulation area into cells and use two arrays HEAD and LIST to store 
linked lists. This method uses no extra memory and has linear performance, i.e., 
0(N) steps. However, the reference locality is broken when molecules in a given cell 
are stored in a long array LIST for a large system. This can be detrimental to the 
performance of the code. 

In our implementation of the MD simulation we build a specially tailored data 



structure to reflect the similarity between the physical problem and the memory 
partition of modern computers. It is shown schematically in Fig. 2. A node, which is 
roughly speaking a piece of computer memory, is created for every particle introduced 
into the system. All the information for that particle is stored in the node. The 
memory of the node is deallocated when the particle moves out of the system. The 
system is divided into cells to reflect the short-ranged nature of the interaction. The 
neighboring environment of every molecule is implemented by a dynamic one-way 
link list in each cell. This one-to-one dynamic management of the simulated systems 
enables a high degree of reference locality and efficient usage of computer resources. 

This data representation also enables layered structures. We already see the 
atom/node layer and the cell layer. Further layers can be built on top of these to take 
into account the structures of e.g., polyatomic molecules or other complex objects. 
Hence this method of constructing the data improves the code reusability and makes 
simulation of complex systems relatively easy. We have, for example, extended our 
code to simulate the flow of a fluid composed of propane (C^Hs) molecules. With 
this data structure, we are able to achieve on a single DECstation 3100 a speed of 
3 x 10~ 5 second/particle/step, comparable or superior to that obtained in somei^ 
but not all0 simulations done on supercomputers for comparable systems. 

III. Results 

All of the figures in this section are for steady-state flow with L = 400cr, w = 
lOOo" = L\, and T = e/k in the source region. The lone exception is the inset in Fig. 5 
which is for a narrow channel with w = 40cj. In all cases the fluid is maintained in 
the source region at a density n = 0.25 /a 2 . The temperature is well above the 
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critical temperature T c ~ 0.55e//c of the two-dimensional Lennard-Jones fluidO and 
the density is well below solid densities n s > 0.8/a 2 so that we will have a one-phase 
system unless sufficient cooling takes place at some point along the channel to produce 
a much lower temperature. Some cooling is expected as the fluid expands and moves 
down the channel; however, collisions with the walls, which are also at T = e/k, will 
compensate somewhat for this effect. 

Figure 3 shows the particle current density along the channel, J x , in units of 
1/(<tt), as a function of x/a. Different sets of points correspond to different intervals 
of y. For example, the points labelled "0-10" represent an average of the current 
density for < y < 10<7, and of the current density for w — Wcr < y < w. Because the 
distribution is symmetric around the line y = w/2, we can average over the data from 
strips symmetrically placed on the two sides of this line. We have done so in order 
to obtain smoother results. The current density averaged across the entire channel is 
also shown. The constancy of this density outside of the source region, as a function 
of x, is an indication that the system is indeed in the steady state. The data in this 
figure, as well as in all of the following ones, are obtained by averaging over a million 
time steps after discarding results from an initial one million steps during which the 
steady state is achieved. Notice also the linear increase of the current densities in the 
source region. 

In Fig. 4 we show the particles' mean velocity in the x direction, u, in units of 
cf/t, as a function of the position across the channel y/w. The velocity is evaluated in 
region II of the channel which is the middle third of the flow area excluding the source 
region. In this part of the channel the end effects are minimized. The velocity u is 
the macroscopic fluid velocity. One may see clearly the velocity slip at the walls. Also 
given in the figure as a solid line is a parabola which has been fit to the data. For an 
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incompressible fluid and laminar flow, the velocity profile is expected to be parabolic. 
From the figure one can see that, even though we have a highly compressible fluid, 
cf. , Figs. 5 and 6, the distribution still fits the parabola quite well. 

The cross-channel density distribution in region II is shown in Fig. 5, and the 
density along the channel, averaged across the channel, appears in Fig. 6. One can 
see that n decreases roughly, but definitely not precisely, linearly as x increases. Also, 
it is is significantly lower in the middle of the channel than toward the sides. The 
cross-channel variation is expected for compressible fluid flow when the velocity is 
large enough.ll! For Mach number less than about 0.3, the density variation is usually 
not great but when it is around 0.5 or larger, which is the case in our system, the 
finite compressibility of the fluid is highly evident. Exact solutions of the hydrody- 
namic equations for the flow of a compressible fluid are known only for a few special 
cases. We propose the following simple arguments to explain qualitatively the density 
distribution observed in our simulations. First, since the local temperature variation 
across the channel is quite small in comparison with the total kinetic energy variation, 
we shall neglect the former. Then to lowest order in the Mach number M we have 
the following generalized Bernoulli equation for isentropic flow, 



where c is the speed of sound; P is the pressure at a point where the flow velocity is 
u; and Po is the pressure where u = 0, or stagnation pressure. The speed of sound 
can be approximated as follows, 



and we can identify AP = P — Pq and Ap = p — po where po is the stagnation density. 



P -P = \ P c 2 M 2 





(3) 
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Then we have, still to second order in M, 

P 1+IM2 i+i^ W 

This function is plotted as a solid line in Fig. 5 with u computed from the simulation. 
For pq we use the density at the downstream end of the source region, 6q = 0.23/ a 2 . 
We did not use the mean density in this region because of the unphysical rescaling 
we use here in the simulations to balance the extra heat generated by bringing in 
new particles. In addition, the speed of sound is temperature-dependent and is not a 
constant throughout the channel leading to a further ambiguity in connection with the 
appropriate choice of c. If we simply choose c to be the ideal gas value c = ^2kT/m 
with T = 0.75/fc which is the average temperature in region II of the channel, cf. , 
fig. 8, that gives the result shown in fig. 5. Qualitatively, the fit is quite reasonable, 
perhaps better than one might have expected. 

Also shown as an inset in Fig. 5 is the cross-channel density distribution in region 
II for a considerably narrower channel with w = 40<r. In this case the density is 
quite constant across the channel except very close to the edges where it decreases 
somewhat. By contrast, the cross-channel velocity distribution for the narrow channel 
shows behavior very similar to that for the side channel except that the peak velocity 
in the center is considerably smaller. The flow is sufficiently slow that no density 
decrease is apparent in the middle of the channel as compared with the density closer 
to the sides at a given x. 

For both the narrow channel and the wide one, there is a significant decrease in 
the particle density very close to the walls, i.e., , within about 2a of the walls. We 
believe this decrease to depend on, and be a consequence of, the manner in which we 
reflect particles from a wall. 
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The local temperature distribution in the channel has been obtained from study 
of the velocity distributions of the particles within each cell. Of course, one may ask 
whether there is a well-defined temperature in an intrinsically non-equilibrium system. 
We obtain the temperature from the equipartition theorem according to which, in 
the absence of any flow, kT = mv x or kT = my 2 where the overlines indicated mean 
values taken over the distribution of particle velocities; both temperatures should be 
the same. However, given that there is some flow velocity u in the x direction as in 
our channel, then if there is a Maxwell velocity distribution centered at this velocity, 
we may introduce temperatures 

kT x = m(v x — u) 2 and kT y = mVy. (5) 

Using these formulas we have computed T x and T y . The cross-channel temperature 
distributions in region II, in units of e/k, are shown in Fig. 7. The two generally agree 
quite well except very close to the edges of the channel where T x is considerably larger 
than T y . One should expect that they would not agree here because the velocity 
distribution is certainly not a moving Maxwell distribution as a consequence of the 
particular manner in which we handle reflections of the particles from the walls. 
Rather, the distribution close to a wall is a superposition of two distributions; the 
reflected particles have a Maxwell distribution which is not moving while the incident 
particles have a net velocity along the x direction. The solid curve in Fig. 7 is the 
best fit of a quartic function to the data. A quartic has been used because that is 
the predicted temperature distribution for an incompressible fluid. One can see that 
the fit is not very good, indicating that the prediction of incompressible fluid theory 
is really not appropriate in this case. 

Figure 8 shows the average temperature distribution along the channel. One 
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can see that T x and T y agree quite well and that the temperature drop, once the 
particles have left the source region, is roughly linear. Several factors contribute to 
the temperature drop. First, the fluid is expanding, which will increase its potential 
energy, and so to the extent that energy is conserved the mean kinetic energy per 
particle must decrease. Second, collisions with the walls offset this effect because 
particles reflected from the walls have on the average a kinetic energy of kT where 
the wall temperature T is the same as the initial temperature of particles injected into 
the source region. Finally, there is also the effect of the increasing flow velocity along 
the channel which has the consequence that, as x increases, more of the kinetic energy 
is going into the flow and correspondingly less energy is available for fluctuations of 
the particles' velocities relative to the flow velocity thereby producing a drop in the 
temperature, as we have defined it, with increasing x. 

IV. Conclusions 

We find that, in two dimensions at least, it is practical to simulate a realistic flow 
system with a source and a drain, and at the same time keep a long enough channel 
that finite size effects associated with the length are not important in regions farthest 
from the source or drain. Sufficiently high Mach numbers to allow for the investigation 
of compressible fluid flow can be achieved as is evidenced by the considerable cross- 
channel density variations observed. The non-uniformity of this density distribution 
can be roughly explained by an equation based on the Bernoulli equation and the 
relation between system compressibility and the sound velocity. This equation also 
predicts, and our results verify, that the density distribution will be quite flat for flow 
with small Mach numbers. 
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We also find that there is good agreement between our results for the flow veloc- 
ity and predictions based on the hydrodynamics of incompressible fluids. Overall the 
results suggest that for fairly high (M ~ 5), but subsonic, Mach number flow, a fairly 
good picture of the system's behavior can be obtained from solutions for incompress- 
ible fluid flow with some corrections of order M 2 to account for the consequences of 
compressibility. 
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Figure Captions 

Fig. 1 Diagram of the two-dimensional channel showing the source and sink regions. 
Walls are indicated by solid lines. Average cross-channel distributions presented 
in some of the following figures are obtained in region II which is equidistant 
from the two ends of the simulation area, excluding the source region. The 
lengths of the channel in the x and y directions are denoted in the text by L 
and w, respectively; the length of the source in the x direction is called L\. 

Fig. 2 Schematic representation of the data structures employed in the molecular dy- 
namics simulations. Each node stores all information such as position, acceler- 
ation, etc., of a particle. The simulation area is divided into rectangular cells. 
Particles (nodes) in each cell are linked by pointers to form a dynamic link list. 

Fig. 3 The particle current density J x , in units of 1/ot and averaged within various 
intervals Ay/w across the channel as indicated by the legends, is shown as a 
function of x/a. Also shown as a solid line is J x averaged across the entire 
channel. These results are for a channel with w = WOa and L = 400<r; the 
source and wall temperature is T = e/k; and the averages are taken over 10 6 
time steps after an initial 10 6 steps during which the system reaches the steady 
state. 

Fig. 4 The steady-state average fluid velocity u in region II is plotted against y/w; u 
is given in units of a/r. The parameters are the same as for Fig. 3. The solid 
line represents the best fit of a parabola to the simulation results. 

Fig. 5 The steady-state cross-channel particle density n in region II in units of 1/a 2 is 
shown as a function of x/a for a wide channel (w = lOOer) as a function of y/w. 
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The solid line is the theory described in the text. The inset displays the particle 
density for a narrow channel (w = 40a). Other parameters are the same as for 
fig. 3. 

Fig. 6 The steady-state average density along the channel in units of I /a 2 is shown as 
a function of x/a. Parameters are the same as in Fig. 3. 

Fig. 7 The steady-state cross-channel temperature distributions T x and T y in region 
II, extracted from the equipartition theorem by computing the particles' kinetic 
energies, are given in units of e/k as functions of y/w. Parameters are the same 
as for Fig. 3. 

Fig. 8 The steady-state temperatures T x and T y in units of e/k, averaged across the 
channel, are shown as functions of x/a. Parameters are the same as for Fig. 3. 
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